\newcommand{\real}{\hbox{\bf R}}


\begin{centering}
\bf Laboratorio de M\'etodos Num\'ericos - Segundo cuatrimestre 2007 \\
\bf Trabajo Pr\'actico N\'umero 3: Homenaje a Los Pumas \\
\end{centering}

\vskip 25pt
\hrule
\vskip 11pt

El objetivo del trabajo pr\'actico es evaluar un m\'etodo para reconstruir
im\'agenes tomogr\'aficas sujetas a ruido, utilizando el m\'etodo de 
aproximaci\'on por cuadrados m\'\i nimos.

\textbf{El m\'etodo de reconstrucci\'on}

El an\'alisis tomogr\'afico de un cuerpo (que suponemos bidimensional y cuadrado) consiste
en emitir se\~nales ac\'usticas que lo atraviesan en diferentes direcciones,
midiendo el tiempo que tarda cada una en atravesarlo. Por ejemplo, la
Figura~1 muestra una posible distribuci\'on de estas se\~nales.

\begin{figure}[h]
	\centering
		\includegraphics[scale=0.4]{tp3.jpg}\\
		\caption{Ejemplo de configuraci\'on de las se\~nales}
	\centering
\end{figure}

\vskip 11pt

Suponemos que el cuerpo est\'a dividido en $n\times n$ celdas cuadradas, de
modo tal que en cada celda las se\~nales ac\'usticas tienen velocidad
constante. Si $d^k_{ij}$ es la distancia que recorre la $k$-\'esima se\~nal
en la celda $ij$ (notar que $d^k_{ij}=0$ si la se\~nal no pasa por esta
celda) y $v_{ij}$ es la velocidad de la se\~nal en esa celda, entonces
el tiempo de recorrido de la se\~nal completa es de
\begin{displaymath}
t_k \ =\ \sum_{i=1}^n \sum_{j=1}^n d^k_{ij} \: v^{-1}_{ij}.
\end{displaymath}
Como resultado del an\'alisis tomogr\'afico, se tienen mediciones del tiempo
que tard\'o cada se\~nal en recorrer el cuerpo. Si se emitieron $m$ se\~nales,
entonces el resultado es un vector $t\in\real^m$, tal que $t_k$ indica el
tiempo de recorrida de la $k$-\'esima se\~nal. Sea $D\in\real^{m\times n^2}$
una matriz cuyas filas se corresponden con las $m$ se\~nales y cuyas columnas
se corresponden con las $n^2$ celdas de la discretizaci\'on del cuerpo, y
tal que la fila $k$ contiene los valores $d^k_{ij}$ en las columnas
correspondientes a cada celda. Entonces, las velocidades de recorrida
originales $v_{ij}$ se pueden reconstruir resolviendo el sistema de
ecuaciones $Ds = t$. El vector soluci\'on $s\in\real^{n^2}$ contiene los
valores inversos de las velocidades originales en cada celda.

Un problema fundamental que debe considerarse es la presencia de errores
de medici\'on en el vector $t$ de tiempos de recorrido. Dado que el sistema
estar\'a en general sobredeterminado, la existencia de estos errores har\'a
que no sea posible encontrar una soluci\'on que satisfaga al mismo tiempo
todas las ecuaciones del sistema. Para manejar este problema, el sistema
$Ds=t$ se resuelve utilizando el m\'etodo de aproximaci\'on por \emph{cuadrados
m\'\i nimos}, obteniendo as\'\i\ una soluci\'on que resuelve de forma aproximada
el sistema original.

\textbf{Enunciado}

El trabajo pr\'actico consiste en implementar un programa que
simule el proceso de tomograf\'\i a y reconstrucci\'on, realizando
los siguientes pasos:
\begin{enumerate}
\item Leer desde un archivo una imagen con los datos reales
(discretizados) de un cuerpo. Para los fines de este procedimiento, se
considera que el valor de cada \emph{pixel} de la imagen corresponde a la
velocidad de recorrida de las se\~nales ac\'usticas al atravesar ese pixel.
\item Ejecutar el proceso tomogr\'afico, generando un conjunto
relevante de se\~nales y sus tiempos de recorrida exactos.
\item Perturbar los tiempos de recorrida con ruido aleatorio.
\item Ejecutar el m\'etodo propuesto en la secci\'on anterior sobre
los datos perturbados, con el objetivo de reconstruir el cuerpo
original. La imagen resultante se debe guardar en un archivo de salida.
\end{enumerate}
El programa debe tomar como par\'ametros los nombres de los archivos
de entrada y de salida, junto con un par\'ametro que permita especificar
el nivel de ruido a introducir en la imagen. El formato de los archivos
de entrada y salida queda a elecci\'on del grupo. Para simplificar la
implementaci\'on, es posible utilizar un formato propio que no sea ninguno
de los formatos est\'andar para guardar im\'agenes\footnote{Un formato para
almacenar im\'agenes que es de f\'acil lectura y tiene cierta difusi\'on es
el formato \textsc{Raw}, ver por ejemplo 
\texttt{http://local.wasp.uwa.edu.au/}$\sim$\texttt{pbourke/dataformats/bitmaps}.}.

Sobre la base de la implementaci\'on, se pide medir la calidad de la imagen
reconstruida en funci\'on del nivel de ruido. Para medir el error de la
imagen resultante se deber\'a utilizar el \emph{error cuadr\'atico medio},
definido como
\begin{displaymath}
\hbox{ecm}(v,\tilde{v})\ =\ \frac{\sum_{i=1}^n \sum_{j=1}^n (v_{ij} - \tilde{v}_{ij})^2}{n^2},
\end{displaymath}
siendo $v_{ij}$ la velocidad de recorrido de la celda $ij$ del cuerpo
original, y $\tilde{v}_{ij}$ la velocidad en la misma celda del cuerpo
reconstruido.

%\vfil \eject

\textbf{Objetivos adicionales}

En forma optativa, se sugiere probar con distintas estrategias geom\'etricas
para generar las se\~nales ac\'usticas (paquetes de se\~nales radiando de
puntos fijos o puntos de inicio m\'oviles, se\~nales partiendo de los cuatro
lados de la imagen o s\'olo de algunos lados, rango de \'angulos de salida
de las se\~nales, etc.) para buscar la estrategia que minimice el error
de la reconstrucci\'on.

Por otra parte, puede ser interesante medir el n\'umero de condici\'on de la
matriz asociada al sistema de ecuaciones normales que resuelve el problema
de cuadrados m\'\i nimos, para analizar la estabilidad de la resoluci\'on.
En caso de que el sistema tenga un n\'umero de condici\'on alto, se pueden
intentar esquemas de regularizaci\'on para mejorar la estabilidad de la
soluci\'on obtenida.

\vskip 15pt

\hrule

\vskip 11pt

Fecha de entrega: Lunes 5 de Noviembre